library(pscl)
library(car)
library(effects)
library(lmtest)
library(stargazer)
library(MASS)
library(multiwayvcov)
library(margins)
library(ggeffects)
library(marginaleffects)
options(scipen = 8)

#laptop
setwd('/Users/joelsievert/Dropbox/Research/Electoral Institutions/DRDR/journal submissions/JPIPE/final submission/Replication Materials')

b <- read.csv('jpipe_tab3.csv')

b$trend <- as.numeric(as.factor(b$congress))

##################
##pension bills##
#################

m1 <- glm.nb(pension ~    total + private_cmte + vet_rank + chamber_terms + maj_party +
             trend  , b)
m2 <- glm.nb(pension ~ total+  senator*private_cmte  + senator*vet_rank  + chamber_terms + maj_party +
               trend   , b)

lrtest(m1, m2)

m3 <- glm.nb(tot_private ~    total + private_cmte + vet_rank + chamber_terms + maj_party + trend  , b)
m4 <- glm.nb(tot_private ~ total +  senator*private_cmte  + senator*vet_rank  + chamber_terms + maj_party + trend   , b)

lrtest.default(m3, m4)


o1 <- coeftest(m1, vcov = cluster.vcov(m1, cluster = b$state_cong))
o2 <- coeftest(m2, vcov = cluster.vcov(m2, cluster = b$state_cong))
o3 <- coeftest(m3, vcov = cluster.vcov(m3, cluster = b$state_cong))
o4 <- coeftest(m4, vcov = cluster.vcov(m4, cluster = b$state_cong))

stargazer(o1, o2, o3, o4, digits = 2)


#########################
##STRUCTURAL BREAK TEST##
#########################

#test structural break, pension#
pooled <- glm.nb(pension ~    total + private_cmte + vet_rank + chamber_terms + maj_party +
                   trend  , b)
sen <- glm.nb(pension ~    total + private_cmte + vet_rank + chamber_terms + maj_party +
                trend  , subset(b, senator == 1))
house<- glm.nb(pension ~    total + private_cmte + vet_rank + chamber_terms + maj_party +
           trend  , subset(b, senator == 0))

pool.ll <- as.numeric(logLik(pooled))
s.ll <- as.numeric(logLik(sen))
h.ll <- as.numeric(logLik(house))

df.test <- (length(coef(sen)) + length(coef(house))) - length(coef(pooled)) 
test.value <- 2* ((s.ll + h.ll) - pool.ll )
1 - pchisq(test.value, df.test)


#test structural break, private#
pooled <- glm.nb(tot_private ~    total + private_cmte + vet_rank + chamber_terms + maj_party +
                   trend  , b)
sen <- glm.nb(tot_private ~    total + private_cmte + vet_rank + chamber_terms + maj_party +
                trend  , subset(b, senator == 1))
house<- glm.nb(tot_private ~    total + private_cmte + vet_rank + chamber_terms + maj_party +
                 trend  , subset(b, senator == 0))

pool.ll <- as.numeric(logLik(pooled))
s.ll <- as.numeric(logLik(sen))
h.ll <- as.numeric(logLik(house))

df.test <- (length(coef(sen)) + length(coef(house))) - length(coef(pooled)) 
test.value <- 2* ((s.ll + h.ll) - pool.ll )
1 - pchisq(test.value, df.test)

################
##predictions##
###############

#########################
#number of pension bills#
#########################

#veteran rank
m2 <- glm.nb(pension ~ total + senator*vet_rank  + senator*private_cmte  + chamber_terms + maj_party 
                 + trend , b)

p1 <- ggpredict(m2, terms = "vet_rank[1:48]", condition = c(senator = 0, total = 28, maj_party = 1, 
                                                             chamber_terms = 1, trend = 15, private_cmte = 0 ), 
                vcov.type = sandwich:vcovCL, vcov.args = c(cong_state), ci.lvl = 0.9 )


p2 <- ggpredict(m2, terms = "vet_rank[1:48]", condition = c(senator = 1, total = 28, maj_party = 1, 
                                                             chamber_terms = 1, trend = 15, private_cmte = 0 ), 
                vcov.type = sandwich:vcovCL, vcov.args = c(cong_state), ci.lvl = 0.9 )


r <- 1:48
j <- 0.25
par(mar = c(5, 4, 3, 2), bg = "white")

plot(p1$x, p1$predicted, ylim = c(0, 16), xlim = c(1,48), type = 'n', xlab='', ylab='', axes= FALSE )

abline(v = seq(1, 47, 2), h = seq(0, 16, 2), lty = 3, col = 'gray')

lines(p1$x , p1$predicted, lwd = 3)
lines(p1$x , p1$conf.low, lwd = 3, lty = 2)
lines(p1$x , p1$conf.high, lwd = 3, lty = 2)

lines(p2$x , p2$predicted, col = 'gray50', lwd = 3)
lines(p2$x , p2$conf.low, lwd = 3, lty = 2, col= "gray50")
lines(p2$x , p2$conf.high, lwd = 3, lty = 2, col = "gray50")


axis(side = 2, at = seq(0, 16, 2), labels = seq(0, 16, 2), las = 2, lwd = 0 ) 
axis(side = 1, at = seq(1, 47, 4), lwd = 0)

mtext("Veteran Rank", side = 1, line = 2.75)
mtext("Predicted Bills", side = 2, line = 2.5)
mtext("Pension Bills", side = 3, line = 0.5, font = 2, cex = 1.15)
legend('bottomleft', border = "n", bty = 'n', lwd = c(3,3), col = c("black", "gray50"), legend = c("Representative", "Senator"))



#veteran rank
m4 <- glm.nb(tot_private ~ total + senator*vet_rank  + senator*private_cmte  + chamber_terms + maj_party 
             + trend , b)

p1 <- ggpredict(m4, terms = "vet_rank[1:48]", condition = c(senator = 0, total = 28, maj_party = 1, 
                                                            chamber_terms = 1, trend = 15, private_cmte = 0 ), 
                vcov.type = sandwich:vcovCL, vcov.args = c(cong_state), ci.lvl = 0.9 )


p2 <- ggpredict(m4, terms = "vet_rank[1:48]", condition = c(senator = 1, total = 28, maj_party = 1, 
                                                            chamber_terms = 1, trend = 15, private_cmte = 0 ), 
                vcov.type = sandwich:vcovCL, vcov.args = c(cong_state), ci.lvl = 0.9 )


r <- 1:48
j <- 0.25
par(mar = c(5, 4, 3, 2), bg = "white")

plot(p1$x, p1$predicted, ylim = c(0, 24), xlim = c(1,48), type = 'n', xlab='', ylab='', axes= FALSE )

abline(v = seq(1, 47, 2), h = seq(0, 24, 2), lty = 3, col = 'gray')

lines(p1$x , p1$predicted, lwd = 3)
lines(p1$x , p1$conf.low, lwd = 3, lty = 2)
lines(p1$x , p1$conf.high, lwd = 3, lty = 2)

lines(p2$x , p2$predicted, col = 'gray50', lwd = 3)
lines(p2$x , p2$conf.low, lwd = 3, lty = 2, col= "gray50")
lines(p2$x , p2$conf.high, lwd = 3, lty = 2, col = "gray50")


axis(side = 2, at = seq(0, 24, 4), labels = seq(0, 24, 4), las = 2, lwd = 0 ) 
axis(side = 1, at = seq(1, 47, 4), lwd = 0)

mtext("Veteran Rank", side = 1, line = 2.75)
mtext("Predicted Bills", side = 2, line = 2.5)
mtext("Private Bills", side = 3, line = 0.5, font = 2, cex = 1.15)
legend('bottomleft', border = "n", bty = 'n', lwd = c(3,3), col = c("black", "gray50"), legend = c("Representative", "Senator"))

###########################
####COMMITTEE MEMBERSHIP###
###########################

#privatee committee
p1 <- ggpredict(m2, terms = "private_cmte[0,1]", condition = c(senator = 0, total = 28, maj_party = 1, trend = 15.5, 
                                                               chamber_terms = 1,  vet_rank = 27), 
                vcov.type = sandwich:vcovCL, vcov.args = c(cong_state), ci.lvl = 0.90 )


p2 <- ggpredict(m2, terms = "private_cmte[0,1]", condition = c(senator = 1, total = 28, maj_party = 1, trend = 15.5, 
                                                               chamber_terms = 1, vet_rank = 27 ), 
                vcov.type = sandwich:vcovCL, vcov.args = c(cong_state), ci.lvl = 0.90 )


j <- 0.025
par(mar = c(5, 4, 3, 2), bg = "white")

plot(p1$x, p1$predicted, ylim = c(0, 12), xlim = c(0, 1), type = 'n', xlab='', ylab='', axes= FALSE )

abline(v = seq(0, 1, 0.2), h = seq(0,12, 2), lty = 3, col = 'gray')

points(c(0.2-j, 0.2 +j), p1$predicted, pch = c(15,16), cex = 1.5, col = "black" )
segments( x0 = c(0.2-j, 0.2 +j), x1 =  c(0.2-j, 0.2 +j), y0 =  p1$conf.low, y1 = p1$conf.high, lwd = 3)

points(c(0.8-j, 0.8 +j), p2$predicted, pch = c(15,16), cex = 1.5, col = "black" )
segments( x0 = c(0.8-j, 0.8 +j), x1 =  c(0.8-j, 0.8 +j), y0 =  p2$conf.low, y1 = p2$conf.high, lwd = 3)

axis(side = 2, at = seq(0, 12, 2), labels = seq(0, 12, 2), las = 2, lwd = 0 ) 
axis(side = 1, at = c(0.2, 0.8), labels = c("Representative", "Senator"), lwd = 0)
mtext("Predicted Bills", side = 2, line = 2.5)
mtext("Pension Bills", side = 3, line = 0.5, font = 2, cex = 1.15)

legend('topright', border = "n", bty = 'n', lwd = c(3,3), pch = c(16, 15), legend = c("Committee Member", "Not Committee Member"))

##############################
#share of private legisaltion#
##############################

#vetaran rank

#privatee committee
p1 <- ggpredict(m4, terms = "private_cmte[0,1]", condition = c(senator = 0, total = 28, maj_party = 1, trend = 15.5, 
                                                               chamber_terms = 1,  vet_rank = 27), 
                vcov.type = sandwich:vcovCL, vcov.args = c(cong_state), ci.lvl = 0.90 )


p2 <- ggpredict(m4, terms = "private_cmte[0,1]", condition = c(senator = 1, total = 28, maj_party = 1, trend = 15.5, 
                                                               chamber_terms = 1, vet_rank = 27 ), 
                vcov.type = sandwich:vcovCL, vcov.args = c(cong_state), ci.lvl = 0.90 )

j <- 0.025
par(mar = c(5, 4, 3, 2), bg = "white")

plot(p1$x, p1$predicted, ylim = c(8, 20), xlim = c(0, 1), type = 'n', xlab='', ylab='', axes= FALSE )

abline(v = seq(0, 1, 0.2), h = seq(8,20, 2), lty = 3, col = 'gray')

points(c(0.2-j, 0.2 +j), p1$predicted, pch = c(15,16), cex = 1.5, col = "black" )
segments( x0 = c(0.2-j, 0.2 +j), x1 =  c(0.2-j, 0.2 +j), y0 =  p1$conf.low, y1 = p1$conf.high, lwd = 3)

points(c(0.8-j, 0.8 +j), p2$predicted, pch = c(15,16), cex = 1.5, col = "black" )
segments( x0 = c(0.8-j, 0.8 +j), x1 =  c(0.8-j, 0.8 +j), y0 =  p2$conf.low, y1 = p2$conf.high, lwd = 3)

axis(side = 2, at = seq(8, 20, 2), labels = seq(8, 20, 2), las = 2, lwd = 0 ) 
axis(side = 1, at = c(0.2, 0.8), labels = c("Representative", "Senator"), lwd = 0)
mtext("Predicted Bills", side = 2, line = 2.5)
mtext("Private Bills", side = 3, line = 0.5, font = 2, cex = 1.15)

legend('topright', border = "n", bty = 'n', lwd = c(3,3), pch = c(16, 15), legend = c("Committee Member", "Not Committee Member"))

